A cluster Monte Carlo algorithm with a conserved order parameter 
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We propose a cluster simulation algorithm for statistical ensembles with fixed order parameter. 
We use the tethered ensemble, which features Helmholtz's effective potential rather than Gibbs's free 
energy, and in which canonical averages are recovered with arbitrary accuracy. For the D = 2, 3 Ising 
model our method's critical slowing down is comparable to that of canonical cluster algorithms. Yet, 
we can do more than merely reproduce canonical values. As an example, we obtain a competitive 
value for the 3D Ising anomalous dimension from the maxima of the effective potential. 
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Monte Carlo simulations fT] constitute one of the most 
important modern tools of theoretical physics. In situa- 
tions that defy an analytical treatment, a Monte Carlo 
computation succeeds by vifandering randomly across the 
system's configuration space. Over the years many meth- 
ods have been proposed in order to optimize this sam- 
pling process. Here we combine tvi^o: fixing some global 
parameter in order to guide the exploration of phase 
space and the use of cluster update algorithms. 

Cluster methods appear in the 1980s '31 as an answer 
to the problem of critical slowing down [3j, the drastic 
deceleration of the dynamics in the neighborhood of the 
critical point. For a system of linear size L, the charac- 
teristic times are t cx (typically 2^2). Only in very 
few situations can one define an efficient cluster method, 
capable of achieving z < 1. 

Beyond the canonical ensemble setting, considering 
global conservation laws is often useful (as in micromag- 
netic or microcanonical Q ensembles). It is known 
that z — 4 — rj for locally conserved order parameter 
dynamics in Ising models [6| (77 is the anomalous dimen- 
sion) . For non-local conservation laws z is smaller 0, S| ■ 
Despite continued research on cluster methods , the de- 
velopment of an efficient cluster method in this situation 
has long been considered somewhat of a challenge Q. 

Here we present a working cluster algorithm with a 
globally conserved order parameter. We employ the 
Tethered Monte Carlo (TMC) framework [io|, which we 
briefly review. TMC is a general approach to reconstruct 
the effective potential [l,[lH. We demonstrate our cluster 
method in the standard benchmark of the D = 2,3 Ising 
model. In the first case it outperforms the Metropolis 
version of , while in the second it exhibits a dynamic 
critical exponent compatible with that of the canonical 
Swendsen-Wang algorithm. 

The tethered ensemble is similar to the micromagnetic 
one, but instead of fixing the magnetization we couple 
it to a Gaussian 'magnetostat' in order to define a new 
parameter, m. This ensemble arises from the canonical 
one through a Legendre transformation that replaces the 
magnetic field h by m. Thus, Helmholtz's effective po- 



tential takes the place of Gibbs's free energy. The main 
observable is the tethered magnetic field h, considered as 
a function of m. 

We shall work on the Ising model in a cubic lattice 
of size N ^ L° and periodic boundary conditions, with 
partition function ((•, •): nearest neighbors) 



Z = ^ exp j3 ^ (Tx<Jy 



(7^ = ±1. 



We shall consider its energy and magnetization, 
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E = Ne 
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We use lowercase for densities so that, for instance, e is 
the energy per bond. Canonical averages are denoted by 
(•)/3, as in the specific heat and susceptibility: 

C^N[{e'),^{e)l], x^N[{m^)i,-{m)l]. (3) 

Note that the probability density (pdf) pi{m) is the 
sum of -I- 1 Dirac deltas. We smooth it by coupling m 
to A^ Gaussian demons to build the tethered ensemble 



M = Nm = M + 
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As the (pi are independent, m ~ m + 1/2. The defini- 
tions of the pdf p(to) for rh and of the effective potential 
f2]\[{rh,P) are straightforward 
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Our use of demons is reminiscent of Creutz's microcanon- 
ical algorithm [13], but we shall integrate the (jji out in 
order to define our tethered averages, 

{O)rn,0 = V- nr~~f — ' 
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where O represents a generic observable and 
u;m{P, rh- {aj) = c-'3^+*^-^^(m - m)(^-2)/2^(^ _ 

(7) 

[9 is Heaviside's step function). Then wc find that 

M-M dm 

Therefore, we can construct the effective potential by in- 
tegrating {h)rh,i3 on rh. Once we have f2N{rh,/3) we can 
compute canonical averages for any given value of the 
external magnetic field h with the formula 



J dm e^t'^"("'^)+''"l(0)A,/3 

/ dm e^[i^N{m,l3)+hrh] 



(9) 



The TMC simulation algorithm consists of four steps: 
(1) Select an appropriate sampling rfii, z = 1, . . . , for 
rh, keeping in mind that m ~ m + 1/2. (2) Run inde- 
pendent simulations for each rhi, measuring the tethered 
averages {0)rh,i3- (3) Integrate (a smooth interpolation 
of) (h)m,i3 to obtain f2N{m,(3). (4) Use equation ([9|) to 
recover the canonical averages. Figure [1] illustrates this 
process for the 3Z3 Ising model. Let us remark that this is 
the algorithm for reproducing canonical averages. How- 
ever, see below, one can also obtain physically relevant 
results directly from the tethered averages. 

In 10] we implemented this algorithm, using Metropo- 
lis dynamics for step (2). Surprisingly enough, magnetic 
observables such as ft. or m presented no critical slowing 
down (other quantities such as the energy presented the 
z « 2 behavior typical of a local update [7|). 




FIG. 1: (color online). Computation of the effective potential 
f2]v(m,/3c), eq. ©, for an L = 128 Ising 3D lattice. Top: We 
simulate each point independently to measure {/i)m,/3^ , eq. ((Sjl , 
(the lines are cubic splines and the errors are much smaller 
than the points). Bottom: This curve is then integrated to 
yield p(m) , which we show in a linear (□) and in a logarithmic 
(O) scale. Inset: (h)rh.,(3^ for m ~ m -|- 1/2 > 1/2. 



Our cluster algorithm, a tethered version of Swendsen- 
Wang, is best explained using the bond-occupation vari- 
ables n^y (=0,1) and the conditional probability dis- 
tributions of [13|. The lattice bond joining neighboring 
sites X and y is occupied if rixy = 1. The occupied 
bonds partition the lattice in connected clusters of size 
Ni, i = 1, . . . ,Nc. Plugging in ^ the identity 



(10) 



where p = 1 — e~ ' , we immediately read the conditional 
probabilities of the spins given the bonds and vice versa: 
(a) As in a canonical case, given the {ctx}, bonds are 
independent and n^y is 1 with probability pS^^^^y- (b) 
Given the n^y the Ni spins within cluster i are equal to 
Si = ±1. The probability of the 2^<^ {Si} configurations 
depends on AI = Y.iti S^N, through eq. jT]): 



pi{S,})<xe'"-^'{M--M) 



\(N-2)/2 



e{M-M). (11) 



Our cluster update will consist, then, in a cluster trac- 
ing using conditioned probability (a) and in a cluster flip- 
ping using (b). For this last step we perform a dynamical 
Monte Carlo, taking the Si ai t = Q from the initial spin 
configuration. At each t we select N'^ <C Nc clusters 
{N'^ K, 5 works just fine). We randomly pick lattice sites, 
selecting the cluster to which they belong, until we find 
N'^ different clusters. We then use (|TT]) to perform a heat 
bath among the 2^'^ configurations with the remaining 
N' clusters fixed 



22l |. We take iVrep such steps 



Taking A'^rep > 1 steps is advantageous because over a 
large number of repetitions many of the Si will eventu- 
ally be flipped, decorrelating the system. Furthermore, it 
takes much longer to trace the clusters than to flip them 
once, so A'lcp can be made relatively large without no- 
ticeably increasing the simulation time. In Figure [2] we 
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FIG. 2: (color online). Spin overlap, eq. (|12p . in a 2D critical 
L = 512 Ising model at the central minimum (□) and right 
maximum (Q) of the effective potential 
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D = 3 
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Cluster 


Met. + Cluster 


L 


Cluster 


Met. + Cluster 


Swendsen-Wang 


16 


2.310(14) 


0.775(3) 


16 


2.135(13) 


0.782(3) 


5.459(3) 


24 


2.440(26) 


0.920(4) 


32 


2.80(3) 


1.134(5) 


7.963(9) 


32 


2.758(20) 


1.055(5) 


48 


3.467(28) 


1.427(8) 


9.831(9) 


64 


3.347(22) 


1.417(7) 


64 


3.88(3) 


1.700(10) 


11.337(12) 


128 


4.11(5) 


1.861(12) 


96 


4.79(5) 


2.152(14) 


13.90(3) 


256 


4.87(4) 


2.391(16) 


128 


5.46(6) 


2.566(17) 


15.90(5) 


512 


5.79(8) 


3.040(24) 


192 


6.54(11) 


3.32(4) 


19.10(9) 


1024 


6.78(8) 


3.70(4) 


256 


7.51(13) 


3.85(5) 


21.83(10) 


Ze 


0.241(7) 




Ze 


0.472(8) 


0.591(4) 


0.460(5) 


xVd.o.f. 


0.36/2 




xVd.o.f. 


5.85/5 


4.61/5 





TABLE I: Integrated autocorrelation times for the energy at m = 0.5 and /3 = /3c for the _D = 2, 3 Ising model. We compare 
the cluster and mixed versions of our TMC algorithm. We also include the results of 14] for canonical Swendsen-Wang. Our 
values for ze are fits to te = AL^'^, where the smallest lattice in range is Lmin = 128 {2D) and Lmin = 32 (3D). 



represent the overlap 

- W^,,])^,, ^^^^ 

which vanishes for completely uncorrelated configura- 
tions. Clearly, the configuration can significantly evolve 
for a fixed distribution of the bonds. A major error re- 
duction (a factor 25 in our largest lattices) is achieved by 
measuring h at each of the N^cp steps. 

Naturally, one must eventually refresh the bond con- 
figuration. A simulation at fixed rh then consists of Nmc 
steps vi^hcrc one traces the clusters and performs A'rop it- 
erations of the random walk in the Si space. We have em- 
pirically found that an A^icp that equilibrates the cluster- 
tracing and cluster-flipping times is close to optimal, and 
very easy to find (one only has to scale A^rcp with N, as 
the tracing of clusters is an 0{N In N) operation). For 
N = 256^ spins, this results in Ny-^p « 5 • 10^. 

The efficiency of a MC method is best assessed through 
the equilibrium autocorrelation function [15] for an ob- 
servable O, Co{t) = ([0(0) - (O)][O(0 - (O)]). The 
slowness of the dynamics can be quantified with the 
integrated autocorrelation time tq- Defining po{t) = 
Co{t)/Co{0), To is just the time integral of po{t) in 
(0, oo). We estimate it with the self-consistent window 
method [l5| for our numerical estimate po of po, 

1 ^ 

To^ - + J2po{t), A = Wto. (13) 
t=i 

We typically use W — but have checked that the re- 
sults are consistent for several Ws. Since tq oc L^^^, one 
is interested in the observable with largest zq {E in our 
case, as is typical for cluster methods). 

We have used te to assess our algorithm's performance 
for the Ising model in two (at (3 = (3c = ln(l + \/2)/2) and 
three dimensions (at (3 = 0.22165459 « /3c [16]). For both 
D's, we find that the te are largest at m = 0.5 (the cen- 
tral minimum of p(?7i)), so we report their values at that 



point in Table HI In 2D we obtain ze — 0.241(7), while 
in iD our dynamic exponent is ze = 0.472(8), compati- 
ble with the Swendsen-Wang value of ze = 0.460(5) [l^ . 
We also include our results with a slightly modified algo- 
rithm, where we take two Metropolis steps each cluster 
step. This mixed algorithm has significantly smaller t'e 
in both dimensions. Paradoxically, in 31?, z'e > ze {z'e is 
unmeasurable in 2D), which probably means that larger 
L would be needed to compute it properly. Since for our 
lattice sizes the mixed algorithm fares better, we shall 
use it hereafter. 

As a proof of TMC's accuracy, we have reproduced 
some of the critical Swendsen-Wang simulations of [l3| 
for the 3D Ising model (Table |TT| . In accordance with 
our Te analysis, the errors with both algorithms are com- 
parable. In 2D, the new cluster algorithm outperforms 
Metropolis. For instance, for an L = 1024 lattice, taking 
10 times fewer steps than in [IG] and the same rh grid, 
cluster errors are 8 times smaller for E, 6 times smaller 
for C and 10 times smaller for x and i^. 

MCS -{e)p C X C ~ 

SW 48 X 10*^ 0.3309822(16) 22.155(18) 21193(13) 82.20(3) 
TMC 50 X 10" 0.3309831(15) 22.174(13) 21202(13) 82.20(6) 

TABLE 11: Comparison of canonical Swendsen-Wang (data 
from 14]) with TMC for an iV = 128^ lattice at /3c. We take 
10® MC steps at each of the 50 points of our m grid (Figure[l}. 
This results in a similar number of MCS for both simulations. 
(5: second-moment correlation length pD. [l7|). 

Yet TMC can do more than reproduce canonical aver- 
ages. Let us compute the anomalous dimension rj. Finite- 
size scaling [l^, lHj tells us that the right maximum of 
p{rh; L) at /3c (Figure [T|) scales as 

mpcak - i = AL-(''+^-2)/2 ^ ^ ^ ^ ^ (A -const.) (14) 

where the dots stand for scaling corrections. 

Unlike in a canonical simulation, we locate rhpeak 
through (ft.)mp„ak,/3 = 0. We simulate two very close val- 
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L\Lmin 


MCS 




V 


xVd.o.f. 


16 


1.0 X 10* 


0.33421(5) 


0.03392(21) 


42.6/6 


32 


1.0 X 10* 


0.23377(4) 


0.03526(30) 


8.79/5 


48 


1.0 X 10* 


0.18956(4) 


0.0360(5) 


4.24/4 


64 


1.0 X 10* 


0.16341(4) 


0.0368(7) 


1.02/3 


96 


1.0 X 10* 


0.13240(4) 


0.0363(12) 


0.78/2 


128 


1.0 X 10* 


0.114083(24) 


0.0373(19) 


0.31/1 


192 


6.0 X lO'^ 


0.09246(4) 






256 


8.2 X 10" 


0.07959(12) 







TABLE III: Position of the peak of p(m; L) for the D = 3 
Ising model at /3c. In each row we report the result of a fit to 
(|14|) (the first column is also the fit's Lmin). 



ues of rh at either side of the peak, and use a linear 
interpolation. Our 3D results are in Table WUl 

We have found that eq. (fTi]) yields remarkably good 
fits for Lmin = 48. Furthermore, increasing results 
in compatible values of rj, with growing errors. Our pre- 
ferred estimate is = 0.0360(7), where we took the cen- 
tral value from Lmin = 48 and the error from Lmin = 64 
to account for systematic effects. This estimate com- 
pares favorably with the best Monte Carlo computation 
known to us, rj — 0.0362(8) fl^ and is compatible with a 
high-temperature expansion value of p 0.03639(15) [l9l | 
(however, both quoted values [l^, [l^l were computed 
with a perfect action, not in the Ising model). 

In summary, we have shown how models with con- 
served order parameter can be efficiently simulated with 
a cluster method. We work in the tethered ensemble 
framework, which allows us to compute the Helmholtz 
effective potential. The method is tested in the D ~ 2,3 
Ising model. For the computation of canonical expecta- 
tion values in large lattices, our cluster algorithm is no 
less efficient than a canonical one in 313 (in 2D the dy- 
namical exponent z « 0.24 is larger than that for the 
canonical algorithm, but still very small). The tethered 
ensemble permits a very efficient computation of quanti- 
ties such as the maxima of the effective potential, which 
would be extremely costly to reproduce in a canonical 
setting. Our estimate for the anomalous dimension of 
the 3D Ising model compares favorably with all previ- 
ous Monte Carlo computations known to us. We plan to 
further_develop this algorithm to study disordered sys- 
tems [20] and the condensation transition 
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